Spontaneous emissions and thermalization of cold bosons in optical lattices 
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We study the thermalization of excitations generated by spontaneous emission events for cold 
bosons in an optical lattice. Computing the dynamics described by the many-body master equa- 
tion, we characterize equilibration timescales in different parameter regimes. For simple observables, 
we find regimes in which the system relaxes rapidly to values in agreement with a thermal distribu- 
tion, and others where thermalization does not occur on typical experimental timescales. Because 
spontaneous emissions lead effectively to a local quantum quench, this behavior is strongly depen- 
dent on the low-energy spectrum of the Hamiltonian, and undergoes a qualitative change at the 
Mott Insulator-superfluid transition point. These results have important implications for the un- 
derstanding of thermalization after localized quenches in isolated quantum g ases, as well as the 
characterization of heating in experiments. 

PACS numbers: 37.10. Jk, 67.85.Hj, 03.75.Lm, 42.50.-p 



Spontaneous emission is a fundamental source of heat- 
ing in optical dipole potentials [TJ [2], and one of the 
key heating sources in current experiments with cold 
atoms in optical lattices [3[ H] . This heating induces non- 
equilibrium dynamics in which thermalization processes 
are expected to play a key role. Typically it is assumed 
that the energy added to the system will be thermalized, 
causing an effective increase in temperature. But does 
that happen? 

This question is a special case of a fundamental prob- 
lem in many-body quantum mechanics: to what extent, 
and under which conditions, will an isolated system un- 
dergo thermalization when perturbed away from equilib- 
rium, in the sense that at long times the system reaches 
a steady state where simple observables take the same 
values as those of a thermal distribution [SHE] • Recently, 
integrable dynamics for strongly interacting cold gases 
confined to move in ID [9] has been observed, allowing 
the study of systems that either do not thermalize [10] 
or undergo generalized thermalization [TTJ [12] . 

In this Letter we investigate these issues by studying 
dynamics induced by spontaneous emission events (in- 
coherent light scattering) for cold bosons in an optical 
lattice [15] , and identify contrasting parameter regimes 
where (i) certain observables relax over short times to 
thermal values, or (ii) the system relaxes on a short 
timescale to states that are clearly non-thermal, even if 
all atoms remain in the lowest Bloch band. In particular, 
the dynamics depends greatly on the low-energy spec- 
trum of the Hamiltonian because spontaneous emissions 
give rise to a local quench, leading to qualitative changes 
at the supernuid-Mott Insulator phase transition. By 
combining time-dependent density matrix renormaliza- 
tion group (t-DMRG) methods [T4HT7] with quantum 
trajectory techniques [T8il20] , we compute the dynamics 
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FIG. 1: (Color online) (a) Absorption and spontaneous emis- 
sion of a lattice photon leads effectively to localization of sin- 
gle atoms. Tunnelling and interactions between atoms then 
redistribute the energy added to the system, (b) Localization 
of an atom in space corresponds initially to a distribution of 
the atom over the whole Brillouin zone (the tails of the quasi- 
momentum distribution are lifted). Subsequent unitary evolu- 
tion leads to a broadened quasi- momentum distribution, i.e., 
the n q =o peak and the tails decrease, while the small quasi- 
momentum components increase (t-DMRG results, U — 2 J, 
N = 48 particles on M = 48 sites, d t = 6, D = 512). 



in the context of real experiments. These results have im- 
portant implications for the characterization of heating 
in current experiments [21]. In fact, the lack of thermal- 
ization of certain excitations may be exploited to enhance 
the realization of fragile many-body states [22 - 25 , lead- 
ing to greater robustness of quantum simulators [26j [27] . 

In Ref. [4 , a many-body master equation was derived 
to describe spontaneous emission processes for bosonic 
atoms. For far-detuned optical fields, we can adiabat- 
ically eliminate the excited atomic levels, obtaining an 
effective equation for ground-state atoms. When the lat- 
tice spacing a is comparable to or greater than the op- 
tical wavelength of scattered photons, a > A, the dy- 
namics of the many-body density operator p is (H = 1), 
p = —i [H, p] +Cip, where the dissipative term describing 
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scattering of laser photons, denoted L\p is 
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and is a multi-band Bose-Hubbard Hamiltonian 
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Here, the 3D band indices are denoted by k, 1, m, n, and 
6^ ^ is a bosonic creation operator for an atom on site 
i in band n. The dissipative dynamics involves scat- 
tering of photons and (in some cases) transitions be- 
tween Bloch bands with the corresponding rates denoted 
Ikimni while the tunnelling in band m between neigh- 
boring sites is denoted J^ n \ on-site interactions are 
denoted [/(M,™,™) anc [ an onsite potential e\ n \ Note 
that this includes the band energy a/ 71 ) as well as (po- 
tentially) an external trapping potential. The parameters 
can be calculated from the microscopic model by expand- 
ing in a basis of Wannier functions [4 , though care must 
be taken to use properly regularized potentials in evalu- 
ating the interaction matrix elements [/(M,™,™) |28l43Q] . 

In deep optical lattices, transition rates for inter-band 
processes coupling neighboring Bloch bands are sup- 
pressed by the square of the Lamb-Dicke parameter, 
n = 27raT/A, where ar is the trap length for the lowest 
band Wannier function. For typical experiments with lat- 
tice depths around V = SE R [with E R = 4tt 2 H 2 /(2mA 2 ), 
where m is the mass of the atom], the suppression is 
n 2 ~ 0.1. In the usual case of a red-detuned optical 
lattice the dominant dissipative processes are thus intra- 
band processes, which return the atoms to their initial 
Bloch band. Processes accessing higher Bloch bands are 
suppressed by a factor of the order rf or greater, and we 
can write an effective two-band master equation describ- 
ing the dynamics of the density operator for atoms in the 
lattice. In what follows, we set 70000 = 7, so that in the 
Lamb-Dicke limit ?/« 1, 71010 = 7oioi = r / 2 7- We also 
use the symbols U = U^°^ and J = J(0). 

Physical interpretation of the master equation - The 
key effect of spontaneous emissions is that the environ- 
ment (in this case the external electromagnetic field) 
measures the position of an atom, effectively localizing 
the atom undergoing a spontaneous emission event on 
a length scale given by the wavelength of the scattered 
photon. Since the wavelength is comparable to the lat- 
tice spacing A ~ a, this quantitatively corresponds to 
localization of the atom on a single lattice site as de- 
picted in Fig. la. Assuming we do not measure the pho- 
ton field, an initially delocalized state of a single atom, 
e.g., an eigenstate of quasimomentum q on the lattice 
|q = 0) oc X^^°^l vac ) 1S decohered into a mixture of 



different possible atom locations. If the atom remains 
in the lowest Bloch band, then p = |q = 0)(q = 0| — » 



W |vac)(vac|b^ , and the atom is spread over the 

entire Brillouin zone in quasi- momentum space. 

Overview of thermalization - As energy is added to 
the system via the dissipative process, it is often pos- 
tulated that Hamiltonian dynamics should be able to 
thermalize the added energy, in the sense that simple 
observables attain values corresponding to a thermal dis- 
tribution with the new mean energy value (E). However, 
the case of inter-band excitations is a clear counterex- 
ample, as the added energy is of the order of band-gap 
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= uj^ — oj(°\ and cannot be thermalized on 
experimental timescales due to a mismatch with the en- 
ergy scales of the lowest Bloch band uo ^> Z7, J. [40 j. This 
is analogous to the collisional stability of doublon pairs 
demonstrated in recent experiments [31]. 

For intra-band processes, it might be expected that 
the system will thermalize, due to the non-integr ability 
of the Bose-Hubbard model except for very small or large 
values of U , similar to what is expected for a quantum 
quench of U/J [32]. However, it is not clear that this 
readily applies to our situation because a spontaneous 
emission event leads to localization of atoms in a local 
quantum quench with excitations that are very low in 
energy. This may result in a lack of thermalization for 
all values of U given the integrable nature of the spec- 
trum at low energies. Below we find that the relaxation 
timescales and equilibrium values strongly depend on the 
interaction characteristics in the lower band. 

Thermalization after a single intra-band spontaneous 
emission event - In order to characterize the thermal- 
ization process, we first consider the situation where the 
system is in the ground state of the single-band Bose- 
Hubbard model \%jjg) at time t = 0, and then undergoes a 
spontaneous emission occurring at site i. In the sense of 
continuous measurement theory [19], the state prepared 
in this way is 
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We consider a weighted ensemble average over the 
sites i with probabilities of spontaneous emission pi oc 
(i/jglnfl^g), and treat a ID system, where we can use t- 
DMRG methods to propagate the state exactly in time. 
However, we expect the mechanisms for thermalization 
that we discuss to also apply in higher dimensions. Note 
that these DMRG results are converged in the matrix 
product state bond dimension D and the truncation of 
the local dimension, d\. 

Fig. lb shows the typical dynamics after a sponta- 
neous emission spreads a particle over the whole Bril- 
louin zone and increases the kinetic energy E kin . The 
interactions between particles transfer some of this in- 
creased kinetic energy to interaction energy, as shown 
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FIG. 2: (Color online) Time evolution of the kinetic energy af- 
ter a single spontaneous emission event (averaged over jumps 
on all possible sites), (a - b) For a superfluid initial state 
(U — 2J), the kinetic energy relaxes to the equilibrium value 
obtained from a Monte-Carlo calculation . For MI states 
(U — 4J), the energy relaxes, but not to E^. The zero value 
of kinetic energy for this plot is the ground state kinetic en- 
ergy Egl n . (c) The difference of the infinite time value of the 
kinetic energy (obtained from an extrapolation of an expo- 
nential fit) to the Monte-Carlo equilibrium energy. For MI 
states with U/J > 3.37, the difference increases rapidly for 
system sizes M = 24, 48, 96 sites (d) The decay rate extracted 
from the exponential fit as a function of U. (t-DMRG results, 
di = 6, D — 256, 512; error bars represent fitting errors and 
statistical errors from QMC). 



explicitly in Fig. 2a for an initial superfluid (SF) state 
with U = 2 J. At t = 0+, E kin is increased by of the 
order of J over the ground state value, and it then re- 
laxes to lower value over a timescale ~ 5/ J in unitary 
time evolution. We obtain an equilibrium value E k ™ 
from path integral Monte-Carlo (QMC) calculations with 
worm-type updates [33] (here in the implementation of 
Ref. [34] - see Ref. [35] for a recent review of the method 
with applications to cold gases) at finite temperature T, 
fitting T to match the value of energy (E) for t > + . It 
is remarkable that this value corresponds to the equilib- 
rium value reached dynamically within statistical errors, 
indicating thermalization of this quantity. In contrast, 
for an initial Mott Insulator (MI, U = 4J) state, E kin 
relaxes on an only slightly longer timescale to an equilib- 
rium value that clearly does not correspond to a thermal 
distribution at the appropriate value of (E). In fact, in 
this parameter regime, thermally induced coherence in 
the MI leads to a E k ™ being close or even below below 
the value of the ground state kinetic energy [36] . 

In Fig. 2c we compare the extrapolated equilibrium ki- 
netic energy, E^ 1 (obtained from an exponential fit) to 
E k ™ for various system sizes and interaction strengths. 
In the SF regime ID Bose-Hubbard model we observe re- 
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FIG. 3: (Color online) Expectation values of the kinetic en- 
ergy of the lowest 1000 eigenstates as a function of the energy 
E in a system with M — 10 sites and N — 10 particles (ex- 
act diagonalization) . The grey line in the upper plot shows 
the distribution of the energy for a canonical ensemble with 
temperature T fitted to give mean energy E. In the SF, the 
eigenvalue expectations are distributed around the Boltzmann 
distribution, but are far from these values in the MI. The 
lower parts show the corresponding occupation probabilities 
for eigenstates after a single spontaneous emission. 



laxation of E^ n to £^ n , while in the MI regime it is far 
from the thermal value. While we might expect for very 
large values of U that thermalization would become diffi- 
cult as the Bose-Hubbard model approaches an integrable 
system with hard-core bosons [10] [37] , the lack of ther- 
malization for values of U / J immediately above the SF- 
MI transition point (when the gap is about A = J/8) is 
very striking [41 . Before performing these calculations, 
we might have expected a crossover behavior, similar to 
that seen in the relaxation rates, as shown in Fig. 2d from 
exponential fits to the long-time behavior of E kin . 

The key to understanding this behavior lies in the 
fact that the spontaneous emissions give rise to a lo- 
cal quantum quench, which only significantly populates 
low-energy eigenstates of the Hamiltonian. Most of the 
amplitude of the resulting wavefunction is in the ground 
state, as shown in Fig. [3] where we plot occupation prob- 
abilities \c a \ 2 and expectation values of the kinetic energy 
(E a \E kin \E a ) in the lowest 1000 energy eigenstates \E a ). 
We see that E kin grows essentially linearly as a function 
of E a , even for U/J ~ 3 near the phase transition, and 
that these values coincide with E k ™ from a Boltzmann 
distribution. Because of the approximate linear depen- 
dence, averages over states with different \c a \ 2 distribu- 
tions will still be approximately equal. Thus, the long 
time average (E kin ) ^ a \c a \ 2 (E a \E kin \E a ) [5] will cor- 
respond to the value for the Boltzmann distribution. As 
soon as we enter the MI phase, between U/J ~ 3 and 
U / J ^ 3.8, there is a qualitative change in the distribu- 
tion of (E a \E km \E a ) values, as depicted in Fig. [ij after 
which we cannot expect to obtain thermal values. In the 
deep MI, the (E a \E kin \E a ) values are far from E k ™, and 
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FIG. 4: (Color online) Quantum trajectory simulation for 
heating which is switched on for a time t = l/J, as illus- 
trated in (a). Heating rates are rates 7 = 0.02 J, 0.04 J, 0.06 J. 
Results are for system of 48 atoms on 48 sites, the standard 
error of the mean is given as shaded area, (b) The increase 
of the total energy during the heating process for a superfluid 
state (U = 2 J), (c) The increase in kinetic energy during 
the heating and the subsequent relaxation. For superfluid 
initial states, the kinetic energy relaxes to the equilibrium 
value (QMC calculations, dashed lines). For a Mott insulat- 
ing initial state, on the same time-scale, the energy does not 
thermalize. This can be seen in (d) where we plot the differ- 
ence between the actual kinetic energy and the equilibrium 
energy (t-DMRG results, D = 256, di = 6, 500 trajectories). 



correspond to excitation of doublon-holon pairs. 

Note that as with thermalization in any closed quan- 
tum system, the behavior depends on the observable 
considered, and sufficiently complicated or non-local ob- 
servables will never thermalize [5 . When we consider 
the quasi-momentum distribution n q in our system with 
open boundary conditions, we find that for large q > 
l/(2a), values agree well with a thermal distribution on 
timescales tJ ~ 5 in the SF for U > 1. However, long 
wavelength modes around n q =o require much longer re- 
laxation timescales. In the MI, the distribution behaves 
qualitatively differently, in that all values of q show small 
discrepancies from the equivalent thermal values, consis- 
tent with what we observed for the kinetic energy. 

Finite-time light scattering and comparison with ther- 
mal distributions - We now consider a specific exper- 
imental setup in which these effects could be observed. 
As depicted in Fig. [4^i, we consider a situation in which 
the background scattering rate is low, and then a mod- 
erate scattering rate is induced for a short time t = 1/ J 
(e.g., via a weak beam with near-resonant light). We 
then switch this off, and observe how the system ther- 
malizes over a timescale of t ~ 5/ J. We compute the 
dynamics of this process by combining t-DMRG meth- 
ods with quantum trajectory techniques, which after a 



stochastic average allow us to determine the many-body 
dynamics from the master equation. 

As shown in Fig. [4)3, the total energy of the sys- 
tem will increase when 7 7^ 0, but will remain con- 
stant once the dissipation is removed. If we denote 
the kinetic energy term of the Hamiltonian for band 

(n) by = -E„,(*J> •7 (n) &i" )t &i" ) . then the time - 
dependence of the total energy from intra-band processes 

in the lowest Bloch band is given by 
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Because the lowest Bloch band for a single particle has 
energies in the range \-2J^°\ +2 the E kin will typ- 
ically be large and negative in the many-body ground 
state, so that this corresponds to a strong increase in the 
energy of the system, as observed in Fig. 4b. 

Correspondingly, the kinetic energy will increase. In 
Fig. 4c, we plot E kin and E k ™ as a function of time. 
As expected from the arguments above, the kinetic en- 
ergy increases much faster than would be expected from 
a thermal distribution with the same increase in total en- 
ergy, and how far away the values are from a thermal dis- 
tribution is dependent on 7. Note that in the experiment 
of Ref. [24], 7 « 0.02 J. In Fig. 4d we plot E kin - E k ™ for 
different values ofU/J. We see clearly that as in the case 
of a single spontaneous emission event above, the kinetic 
energy will relax towards the expected equilibrium values 
in the superfluid regime. Strikingly, this is not the case in 
the Mott Insulator, where the system remains well away 
from the equilibrium value on the timescales calculated. 
Note that while here the energy increase is small, as we 
use parameters where few spontaneous emission events 
occur to allow quantitative numerical treatments, exper- 
iments could work with larger increases through faster 
scattering rates or longer excitation timescales. 

Conclusions - We showed that for bosons in an op- 
tical lattice, a qualitative change in the thermalization 
behavior after spontaneous emission events will occur at 
the SF-MI transition point. Certain simple quantities 
including the kinetic energy and quasi-momentum distri- 
bution settle rapidly to a steady state. However, while in 
some cases these values correspond to a thermal distribu- 
tion, in others the values are demonstrably non-thermal. 
These findings, presented here for a uniform system, re- 
main valid in the presence of a harmonic trap, as is shown 
by results presented in the supplementary material. The 
lack of complete thermalization implies that the specific 
effects on specific many-body states must be considered. 

For the realization of low-temperature states in optical 
lattices, this may lead in some regimes to a greater ro- 
bustness of the states being produced, as thermalization 
of all the energy added in a spontaneous emission event 
would often lead to temperatures above those required for 
realization of fragile types of order [22 -25 . Instead, be- 
cause the dynamics must be treated as a non-equilibrium 
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situation on a case-by-case basis, much of the interesting 
order in the state can survive on substantial timescales. 

We thank I. Bloch, D. Boyanovsky, W. Ketterle, 
S. Langer, H. Pichler, U. Schneider, D. Weiss, and 
P. Zoller for helpful and motivating discussions. This 
work was supported in part by AFOSR grant FA9550-12- 
1-0057, and by a grant from the US Army Research Office 
with funding from the DARPA OLE program. Compu- 
tational resources were provided by the Center for Sim- 
ulation and Modeling at the University of Pittsburgh. 



[1] J. P. Gordon and A. Ashkin, Phys. Rev. A 21, 1606 
(1980). 

[2] J. Dalibard and C. Cohen-Tannoudji, Journal of Physics 

B: Atomic and Molecular Physics 18, 1661 (1985). 
[3] F. Gerbier and Y. Castin, Phys. Rev. A 82, 013615 

(2010) . 

[4] H. Pichler, A. J. Daley, and P. Zoller, Phys. Rev. A 82, 
063605 (2010). 

[5] M. Rigol, V. Dunjko, and M. Olshanii, Nature 452, 854 
(2008). 

[6] M. A. Cazalilla, R. Citro, T. Giamarchi, E. Orignac, and 
M. Rigol, Rev. Mod. Phys. 83, 1405 (2011). 

[7] M. Rigol and M. Srednicki, Phys. Rev. Lett. 108, 110601 
(2012). 

[8] M. Srednicki, Phys. Rev. E 50, 888 (1994). 
[9] T. Kinoshita, T. Wenger, and D. S. Weiss, Nature 440, 
900 (2006). 

[10] M. Rigol, Phys. Rev. Lett. 103, 100403 (2009). 

[11] A. C. Cassidy, C. W. Clark, and M. Rigol, Phys. Rev. 

Lett. 106, 140405 (2011). 
[12] M. Rigol and M. Fitzpatrick, Phys. Rev. A 84, 033640 

(2011) . 

[13] I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 
80, 885 (2008). 

[14] G. Vidal, Phys. Rev. Lett. 93, 040502 (2004). 

[15] A. J. Daley, C. Kollath, U. Schollwock, and G. Vidal, 
Journal of Statistical Mechanics: Theory and Experi- 
ment p. P04005 (2004). 

[16] S. R. White and A. E. Feiguin, Phys. Rev. Lett. 93, 
076401 (2004). 

[17] F. Verstraete, V. Murg, and J. I. Cirac, Advances in 

Physics 57, 143 (2008). 
[18] K. M0lmer, Y. Castin, and J. Dalibard, J. Opt. Soc. Am. 

B 10, 524 (1993). 
[19] C. W. Gardiner and P. Zoller, Quantum Noise (Springer, 

Berlin, 2005). 

[20] H. J. Carmichael, An Open Systems Approach to Quan- 
tum Optics (Springer, Berlin, 1993). 

[21] D. C. McKay and B. DeMarco, Reports on Progress in 
Physics 74, 054401 (2011). 

[22] S. Fuchs, E. Gull, L. Pollet, E. Burovski, E. Kozik, T. Pr- 
uschke, and M. Troyer, Phys. Rev. Lett. 106, 030401 
(2011). 

[23] R. Jordens, L. Tarruell, D. Greif, T. Uehlinger, 
N. Strohmaier, H. Moritz, T. Esslinger, L. De Leo, 
C. Kollath, A. Georges, et al., Phys. Rev. Lett. 104, 
180401 (2010). 

[24] S. Trotzky, L. Pollet, F. Gerbier, U. Schnorrberger, 



I. Bloch, N. V. Prokof'ev, B. Svistunov, and M. Troyer, 

Nat Phys 6, 998 (2010). 
[25] T. Esslinger, Annual Review of Condensed Matter 

Physics 1, 129 (2010). 
[26] I. Bloch, J. Dalibard, and S. Nascimbene, Nat Phys 8, 

267 (2012). 

[27] J. I. Cirac and P. Zoller, Nat Phys 8, 264 (2012). 
[28] H. P. Biichler, Phys. Rev. Lett. 104, 090402 (2010). 
[29] M. J. Mark, E. Haller, K. Lauber, J. G. Danzl, A. J. 
Daley, and H.-C. Nagerl, Phys. Rev. Lett. 107, 175301 

(2011) . 

[30] M. J. Mark, E. Haller, K. Lauber, J. G. Danzl, 
A. Janisch, H. P. Biichler, A. J. Daley, and H.-C. Nagerl, 
Phys. Rev. Lett. 108, 215302 (2012). 

[31] N. Strohmaier, D. Greif, R. Jordens, L. Tarruell, 
H. Moritz, T. Esslinger, R. Sensarma, D. Pekker, E. Alt- 
man, and E. Dernier, Phys. Rev. Lett. 104, 080401 
(2010). 

[32] G. Biroli, C. Kollath, and A. M. Lauchli, Phys. Rev. Lett. 

105, 250401 (2010). 
[33] N. Prokof'ev, B. Svistunov, and I. Tupitsyn, Journal of 

Experimental and Theoretical Physics 87, 310 (1998), 

ISSN 1063-7761. 
[34] L. Pollet, K. V. Houcke, and S. M. Rombouts, Journal 

of Computational Physics 225, 2249 (2007), ISSN 0021- 

9991. 

[35] L. Pollet, Reports on Progress in Physics 75, 094501 

(2012) . 

[36] E. Toth and P. B. Blakie, Physical Review A 83, 
021601(R) (2011). 

[37] M. Rigol, V. Dunjko, V. Yurovsky, and M. Olshanii, 
Phys. Rev. Lett. 98, 050405 (2007). 

[38] T. Stoferle, H. Moritz, C. Schori, M. Kohl, and 
T. Esslinger, Phys. Rev. Lett. 92, 130403 (2004). 

[39] P. Barmettler, D. Poletti, M. Cheneau, and C. Kollath, 
Physical Review A 85, 053625 (2012). 

[40] Note that collisional processes between two or more 
atoms in the first excited band can return particles to 
the lowest band while exciting atoms to higher bands. 
This doesn't affect the conclusion that the bandgap en- 
ergy cannot be thermalized with the atoms in the lowest 
band. 

[41] Although from our calculations we cannot rule out a sec- 
ond relaxation process to a thermal distribution for much 
larger systems or on much longer timescales, it is clear 
that a striking qualitative change in behavior occurs here, 
leading to a lack of thermalization on typical experimen- 
tal timescales 



6 



SUPPLEMENTAL MATERIAL 
EXCITATIONS TO HIGHER BANDS 



The intra-band processes described in the main text will dominate the dynamics for red-detuned optical lattices, 
and spontaneous emissions induced by near-resonant light with a constant intensity profile. This dominance could be 
increased by adding an additional step in which the lattice is made very deep, and the interactions very weak during 
the process of scattering near-resonance light. This decreases rf , and thus the probability of inter-band processes. 
Band-mapping techniques are also an important diagnostic tool to determine the occupation of higher bands [38] . 

Here we consider also the case of a spontaneous emission event, which takes an atom to the first excited band and 
study the time evolution of the state 

under the two-band Bose-Hubbard model, i.e., the Hamiltonian H = Hq + Hi + i?oi? where i^o,i is the respective 
single band Bose-Hubbard Hamiltonian for each band, and H Q1 = (C/oi/2) ^ with Uqi = JJ 1,0,1,0 . 

Since we consider a transition to the first excited band, the tunneling amplitude in this band < 0, and we consider 
a shallow lattice with Vq = 5i^R, so that ~ — 7 J. Furthermore, typically Uoi ~ U and we choose the two to be 
equal. After the spontaneous emission (Hi) = 0, since the particle in the excited band is localized and there is no 
interaction. As shown in Fig. [5^i, we see that in the subsequent evolution the particle in the excited band starts to 
delocalize, thereby lowering its energy. Through collisions with particles in the lowest band, energy is transferred to 
the lowest band and (Hq) increases. On the timescale of our simulation no steady state is achieved. 



EFFECTS OF A HARMONIC TRAP 



In a realistic experimental setup, the particles will always be confined by a harmonic trap. This can lead to 
situations, in which parts of the system have superfluid components despite the fact that U/J> 3.37. In Fig. ^jp we 
show results for the evolution after a single spontaneous emission events, i.e., of the state Q in the presence of a 
harmonic confinement = ojt — ^o) 2 (center site zq = (M + l)/2) with ujt = 0.012 J for a system with M = 48 
sites. We again average over jumps on all possible sites weighted with the probabilities pi oc (^|n 2 |^). In the upper 
panel we find thermalization for the kinetic energy in the superfluid state with U = 2 J. As in the case of box boundary 
conditions, the kinetic energy relaxes to a value corresponding to a thermal ditribution on the experimentally relevant 
time-scales. For larger interactions (U = 5J), we see from the density profile (insets) that the system is not in a MI 
state with unit filling (these appear only for values around U/J ~ 10). Nevertheless, in this case we find that the 
system does not relax to a value of a thermal distribution. Note that the Monte-carlo thermal kinetic energy in this 
case is in fact above the value after the jump. In the case of a Mott insulator in a trap (U = 10 J) with unit filling 
(seen by the density profile in the inset), we find that the system does not approach a steady state kinetic energy but 
shows oscillations, which can be explained by boundary effects of deflected doublon-hole pairs (see below). 



DOUBLON HOLON RESULTS FOR STRONG INTERACTIONS 



In the limit of strong interactions we can show that there cannot be any thermalization even far away from the regime 
where the system assumes an ideal or a hard-core boson Mott-insulator [TO] [37] . Therefore we perform calculations 
in a regime where we can restrict the occupations of each site to 0,1, and 2 bosons. In this regime one can map 
the system to a model of free fermionic doublon-holon pairs. Using a generalized Jordan- Wigner transformation as 
described in [39] , it is straightforward to show that the time-evolved state after the spontaneous emission on site m is 

x e-^^^^^^c^^c^Jvac). (7) 

Here, the state |vac) is the vacuum, which is annihilated by the free quasi-particles Cd and c^, which are superpositions 
of double occupated sites and holes. w(q,q') is a characteristic wave-packet (example shown in Fig(5b), and the 
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FIG. 5: (Color online) (a) Time evolution of the energy in the two lowest bands after a single spontaneous emission, which 

2)(b) 



[i] 



j[2] 



takes a particle to the excited band (averaged over all sites) (M = N — 24, t-DMRG results, D — 256, d\ 
Relaxation in the lowest band in the presence of a harmonic trap (u)j> — 0.012 J, N — 48, M — 64). The insets show the initial 
density profiles of the ground-states in the trap. (t-DMRG results, D = 512, di = 6, 8) (c) The doublon-holon wave packet 
created by a single spontaneous emission in the center of the system for U = 10J. (d) Time-dependent double occupancy, 
which reveals the propagation of the wave-packet (dashed white line: analytical result) (e) Results for the time-evolution of the 
kinetic energy for strong interactions. The inset shows the evolution after the spontaneous emission. Boundary effects become 
important once the wave-packet reaches the boundary. The kinetic energy is nearly constant as function of time and matches 
the analytical result (grey line) (t-DMRG results, D = 512, di = 6). 



dispersion relations for c^/h are denoted by td,h(q)- As shown in Fig. 5d, after a spontaneous emission the free 
particles spread through the system with a characteristic group velocity before boundary effects become important. 

In addition, from state ^ it is straightforward to show, that any observable in quasi- momentum space, i.e., any 
operator, which can be written in the form O = O q such as the kinetic energy is time-independent, and therefore 
no relaxation is possible in this regime. In Fig. [5^ we compare the analytical result for the difference of the kinetic 
energy to the ground-state to a numerical t-DMRG simulation after a single spontaneous emission in the center of the 
system. We find that in agreement to the doublon-holon calculation, the kinetic energy after the jump remains nearly 
constant (up to some small initial increase due to higher order effects) as long as boundary effects become important 
at t J ~ 6. , 



